/*
 * 	This file is a modification of phylip's drawtree.c. It is also part 
 *  of MVD_GUI.
 *
 *  MVD_GUI is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MVD_GUI is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with MVD_GUI.  If not, see <http://www.gnu.org/licenses/>.
 *  This version copyright 2010 Desmond Schmidt
 */
/* In addition to the above, the following notice also applies */
/* Version 3.6.  Copyright (c) 1986-2004 by the University of Washington and
  Written by Joseph Felsenstein and Christopher A. Meacham.  Additional code
  written by Sean Lamont, Andrew Keefe, Hisashi Horino, Akiko Fuseki, Doug
  Buxton and Michal Palczewski.  Permission is granted to copy, distribute,
  and modify this program provided that (1) this copyright message is
  not removed and (2) no fee is charged for this program. */

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <syslog.h>
#include <stdlib.h>
#include "phylip.h"
#include "draw.h"

#define GAP 0.5
#define MAXITERATIONS 100
#define MINIMUMCHANGE   0.0001
#ifndef INFINITY
#define INFINITY 0x7fffffff
#endif


FILE *plotfile;
FILE *intree;
char plotFilePath[FILE_NAME_LEN];
char inTreePath[FILE_NAME_LEN];
int maxNumOfIter;
int numlines;
int nextnode;
node *root;
node *grbg;
node *where;
pointarray nodep;
boolean goteof;
boolean haslengths;
boolean rotate;
boolean improve;
boolean nbody;
boolean regular;
boolean rescaled;
boolean labelavoid;
boolean uselengths;
double topoflabels;
double bottomoflabels;
double rightoflabels;
double leftoflabels;
double maxx;
double maxy;
double minx;
double miny;
double xmargin;
double xsize;
double bscale;
double ymargin;
double ysize;
double maxchange;
double xcorner;
double ycorner;
double xunitspercm;
double yunitspercm;
double labelheight;
double charht;
double yscale;
double expand;
double pagex;
double pagey;
double paperx;
double papery;
double hpmargin;
double vpmargin;
double xscale;
double xoffset;
double yoffset;
double treeangle;
double ark;
double labelrotation;
double *textlength, *firstlet;
plottertype plotter;
labelorient labeldirec;
growth grows;
char fontname[FONT_NAME_LEN];
enum { yes, no } penchange, oldpenchange;
struct stackElem
{
  /* This is actually equivalent to a reversed link list; pStackElemBack
     point toward the direction of the bottom of the stack */
  struct stackElem *pStackElemBack;
  node *pNode;
};
typedef struct stackElem stackElemType;
/* calculate force acting between 2 nodes and return the force in pForce. 
   Remember to pass the index subnodes to this function if needed. 
   Force should always be positive for repelling.  Angle changes to 
   indicate the direction of the force.  The value of INFINITY is the cap
   to the value of Force.  
   There might have problem (error msg.) if pFromSubNode and pToSubNode 
   are the same node or the coordinates are identical even with double 
   precision.  */
void force_1to1(node *pFromSubNode, node *pToSubNode, double *pForce,
	double *pAngle, double medianDistance)
{
  	double distanceX, distanceY, distance, norminalDistance;
  	distanceX = pFromSubNode->xcoord - pToSubNode->xcoord;
  	distanceY = pFromSubNode->ycoord - pToSubNode->ycoord;
  	distance = sqrt( distanceX*distanceX + distanceY*distanceY );
  	norminalDistance = distance/medianDistance;

  	if (norminalDistance < epsilon)
  	{
    	*pForce = INFINITY;
  	}
  	else
  	{
    	*pForce = (double)1 / (norminalDistance * norminalDistance);
    	if (*pForce > INFINITY) *pForce = INFINITY;
  	}
  	*pAngle = computeAngle(pFromSubNode->xcoord, pFromSubNode->ycoord,
    	pToSubNode->xcoord, pToSubNode->ycoord);
  return;
}
/* pToSubNode is where all the relevent nodes apply forces to.
   All branches are visited except the branch contains pToSubNode.
   pToSubNode must be one of the branch out of the current Node (= out of one
   of the subnode in the current subnodes set.)
   Most likely pPivotSubNode is not the index subNode!  In any case, only
   the leafs are consider in repelling force; so, no worry about index subNode.
   pTotalForce and pAngle must be set to 0 before calling this function for the
   first time, or the result will be invalid.  pPivotSubNode is named for 
   external interface.  When calling totalForceOnNode() recursively, 
   pPivotSubNode should be thought of as pFromSubNode.
 */
void totalForceOnNode(node *pPivotSubNode, node *pToSubNode,
	double *pTotalForce, double *pAngle, double medianDistance)
{
  	node *pSubNode;
  	double force, angle, forceX, forceY, prevForceX, prevForceY;

  	pSubNode = pPivotSubNode;

  	/* visit the rest of the branches of current node; the branch attaches to 
     	the current subNode may be visited in the code down below. */
  	while (pSubNode->next != NULL && pSubNode->next != pPivotSubNode)
  	{
    	pSubNode = pSubNode->next;
    	if ( pSubNode->back != NULL && pSubNode->back != pToSubNode) 
      		totalForceOnNode(pSubNode->back, pToSubNode, pTotalForce, pAngle, 
          	medianDistance);
  	}

  	/* visit this branch; You need to visit it for the first time - at root only!
   	 *
     * Modified so that all nodes are visited and calculated forces, instead of 
     * just the leafs only. 
     * use pPivotSubNode instead of pSubNode here because pSubNode stop short
     *  just before pPivotSubNode (the entry node) */
  	if ( pPivotSubNode == root && pPivotSubNode->back != NULL 
       	&& pPivotSubNode->back != pToSubNode) 
    totalForceOnNode(pPivotSubNode->back, pToSubNode, pTotalForce, pAngle, 
        medianDistance);
  	/* Break down the previous sum of forces to components form */
  	prevForceX = *pTotalForce * cos(*pAngle);
  	prevForceY = *pTotalForce * sin(*pAngle);
  	force_1to1(nodep[pPivotSubNode->index-1], pToSubNode, &force, &angle, 
      medianDistance);
  	/* force between 2 nodes */
  	forceX = force * cos(angle);
  	forceY = force * sin(angle);
  	/* Combined force */
  	forceX = forceX + prevForceX;
  	forceY = forceY + prevForceY;
  	/* Write to output parameters */
  	*pTotalForce = sqrt( forceX*forceX + forceY*forceY );
  	*pAngle = computeAngle((double)0, (double)0, forceX, forceY);
  	return;
}
/* Return the equivalent value of angle that is within
   0 to 2*pi */
double capedAngle(double angle)
{
  	while (angle < 0 || angle >= 2*pi)
  	{
    	if(angle < 0)
    	{
      		angle = angle + 2*pi;
    	}
    	else if (angle >= 2*pi)
    	{
      		angle = angle - 2*pi;
    	}
  	}
  	return angle;
}
/* it return the sign of the moment caused by the force, applied
   to the tip of the refereceVector; the root of the refereceVector
   is the pivot. */
double signOfMoment(double xReferenceVector, double yReferenceVector,
                double xForce, double yForce)
{
  	double angleReference, angleForce, sign;

  	angleReference = computeAngle((double)0, (double)0, xReferenceVector, 
                                yReferenceVector);
  	angleForce = computeAngle((double)0, (double)0, xForce, yForce);
  	angleForce = capedAngle(angleForce);
  	angleReference = capedAngle(angleReference);

  	/* reduce angleReference to 0 */
  	angleForce = angleForce - angleReference;
  	angleForce = capedAngle(angleForce);

  	if (angleForce > 0 && angleForce < pi)
  	{
    	/* positive sign - force pointing toward the left of the reference 
       	line/vector.  */
    	sign = 1;
  	}
  	else
  	{
    	/* negative sign */
    	sign = -1;
  	}
  	return sign;
}
/* added - danieyek 990205 */
/* pStackTop must be the current top element of the stack, where we pop an
   element from the top of it.  
   ppStackTop must be the location where we can find pStackTop.
   This function "returns" the revised top (element) of the stack through 
   the output parameter, ppStackTop. 
   The last element on the stack has the "back" (pStackElemBack) pointer
   set to NULL.  So, when the last element is poped, ppStackTop will be
   automatically set to NULL.  If popNodeFromStack() is called with 
   ppStackTop = NULL, we assume that it is the error caused by over popping
   the stack.
*/
void popNodeFromStack(stackElemType **ppStackTop, node **ppNode)
{
  	stackElemType *pStackT;

  	if (ppStackTop == NULL)
  	{
    	syslog(LOG_ERR,"drawtree - a call to pop while the stack is empty.\n");
    	exit(1);
  	}
  	pStackT = *ppStackTop;
  	*ppStackTop = pStackT->pStackElemBack;
  	*ppNode  = pStackT->pNode;
  	free(pStackT);
  	return;
}
void polartrav(node *p, double xx, double yy, double firstx,
    double firsty, double *leftx, double *lefty, double *rightx,
	double *righty)
{
  	/* go through subtree getting left and right vectors */
  	double x, y, xxx, yyy, labangle = 0;
  	boolean lookatit;
  	node *pp;

  	lookatit = true;
  	if (!p->tip)
    	lookatit = (p->next->next != p || p->index != root->index);
  	if (lookatit) 
	{
    	x = nodep[p->index - 1]->xcoord;
    	y = nodep[p->index - 1]->ycoord;
    	if (p->tip) 
		{
      		if (labeldirec == fixed) 
			{
        		labangle = pi * labelrotation / 180.0;
        		if (cos(p->oldtheta) < 0.0)
          			labangle = labangle - pi;
      		}
      		if (labeldirec == radial)
        		labangle = p->theta;
      		else if (labeldirec == along)
        		labangle = p->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		xxx = x;
      		yyy = y;
      		if (labelavoid) 
			{
        		if (labeldirec == middle) 
				{
          			xxx += GAP * labelheight * cos(p->oldtheta);
          			yyy += GAP * labelheight * sin(p->oldtheta);
          			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			if (textlength[p->index - 1] * sin(p->oldtheta) < 1.0)
            			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			else
            			xxx += 0.5 * labelheight * cos(labangle)
                     		* textlength[p->index - 1];
          			yyy += labelheight * sin(labangle) * textlength[p->index - 1];
        		}
        		else 
				{
          			xxx += GAP * labelheight * cos(p->oldtheta);
          			yyy += GAP * labelheight * sin(p->oldtheta);
          			xxx -= labelheight * cos(labangle) * 0.5 * firstlet[p->index - 1];
          			yyy -= labelheight * sin(labangle) * 0.5 * firstlet[p->index - 1];
          			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			yyy += labelheight * sin(labangle) * textlength[p->index - 1];
        		}
      		}
      		if ((yyy - yy) * firstx - (xxx - xx) * firsty < 0.0) 
			{
        		if ((yyy - yy) * (*rightx) - (xxx - xx) * (*righty) < 0.0) 
				{
          			(*rightx) = xxx - xx;
          			(*righty) = yyy - yy;
        		}
      		}
      		if ((yyy - yy) * firstx - (xxx - xx) * firsty > 0.0) 
			{
        		if ((yyy - yy) * (*leftx) - (xxx - xx) * (*lefty) > 0.0) 
				{
          			(*leftx) = xxx - xx;
          			(*lefty) = yyy - yy;
        		}
      		}
    	}
    	if ((y - yy) * firstx - (x - xx) * firsty < 0.0) 
		{
      		if ((y - yy) * (*rightx) - (x - xx) * (*righty) < 0.0) 
			{
        		(*rightx) = x - xx;
        		(*righty) = y - yy;
      		}
    	}
    	if ((y - yy) * firstx - (x - xx) * firsty > 0.0) 
		{
      		if ((y - yy) * (*leftx) - (x - xx) * (*lefty) > 0.0) 
			{
        		(*leftx) = x - xx;
        		(*lefty) = y - yy;
      		}
    	}
  	}
  	if (p->tip)
    	return;
  	pp = p->next;
  	while (pp != p) 
	{
    	if (pp != NULL)
      		polartrav(pp->back,xx,yy,firstx,firsty,leftx,lefty,rightx,righty);
    	pp = pp->next;
  	}
}
void polarize(node *p, double *xx, double *yy)
{
  	double TEMP, TEMP1;

  	if (fabs(p->xcoord - (*xx)) > epsilon)
    	p->oldtheta = atan((p->ycoord - (*yy)) / (p->xcoord - (*xx)));
  	else if (p->ycoord - (*yy) > epsilon)
      	p->oldtheta = pi / 2;
  	if (p->xcoord - (*xx) < -epsilon)
    	p->oldtheta += pi;
  	if (fabs(p->xcoord - root->xcoord) > epsilon)
    	p->theta = atan((p->ycoord - root->ycoord) / (p->xcoord - root->xcoord));
  	else if (p->ycoord - root->ycoord > 0.0)
      	p->theta = pi / 2;
    else
      	p->theta = 1.5 * pi;
  	if (p->xcoord - root->xcoord < -epsilon)
    	p->theta += pi;
  	TEMP = p->xcoord - root->xcoord;
  	TEMP1 = p->ycoord - root->ycoord;
  	p->r = sqrt(TEMP * TEMP + TEMP1 * TEMP1);
}
/* added - danieyek 990128 */
/* After calling tilttrav(), if you don't polarize all the nodes on the branch
   to convert the x-y coordinates to theta and radius, you won't get result on 
   the plot!  This function takes a subnode and branch out of all other subnode
   except the starting subnode (where the parent is), thus converting the x-y 
   to polar coordinates for the branch only.  xx and yy are purely "inherited"
   features of polarize().  They should have been passed as values not 
   addresses. */
void polarizeABranch(node *pStartingSubNode, double *xx, double *yy)
{
  	node *pSubNode;

  	pSubNode = pStartingSubNode;
  	/* convert the current node (note: not subnode) to polar coordinates. */
  	polarize( nodep[pStartingSubNode->index - 1], xx, yy);

  	/* visit the rest of the branches of current node */
  	while (pSubNode->next != NULL && pSubNode->next != pStartingSubNode)
  	{
    	pSubNode = pSubNode->next;
    	if ( pSubNode->tip != true )
    	  	polarizeABranch(pSubNode->back, xx, yy);
  	}
  	return;
}
/* traverse to move successive nodes */
void tilttrav(node *q, double *xx, double *yy, double *sinphi,
	double *cosphi)
{
  	double x, y;
  	node *pp;

  	pp = nodep[q->index - 1];
  	x = pp->xcoord;
  	y = pp->ycoord;
  	pp->xcoord = (*xx) + (x - (*xx)) * (*cosphi) + ((*yy) - y) * (*sinphi);
  	pp->ycoord = (*yy) + (x - (*xx)) * (*sinphi) + (y - (*yy)) * (*cosphi);
  	if (q->tip)
    	return;
  	pp = q->next;
  	while (pp != q) 
	{
	    /*
	    if (pp != root)
	    */
	    if (pp->back != NULL)
 	     	tilttrav(pp->back,xx,yy,sinphi,cosphi);
    	pp = pp->next;
  	}
}
/* added - danieyek 990204 */
/* pStackTop must be the current top element of the stack, where we add another
   element on top of it.  
   ppStackTop must be the location where we can find pStackTop.
   This function "returns" the revised top (element) of the stack through 
   the output parameter, ppStackTop. 
   The last element on the stack has the "back" (pStackElemBack) pointer
   set to NULL.  So, when the last element is poped, ppStackTop will be
   automatically set to NULL.  If popNodeFromStack() is called with 
   ppStackTop = NULL, we assume that it is the error caused by over popping
   the stack.
*/
void pushNodeToStack(stackElemType **ppStackTop, node *pNode)
{
  	stackElemType *pStackElem;

  	if (ppStackTop == NULL)
  	{
    	/* NULL can be stored in the location, but the location itself can't 
       	be NULL! */
    	syslog(LOG_ERR,"drawtree - error using pushNodeToStack(); " 
        "ppStackTop is NULL.\n");
    	exit(-1);
  	}
  	pStackElem = (stackElemType*)Malloc( sizeof(stackElemType) );
  	pStackElem->pStackElemBack = *ppStackTop;
  	/* push an element onto the stack */
  	pStackElem->pNode = pNode;
  	*ppStackTop = pStackElem;
  	return;
}
/* polar coordinates of a node relative to start */
void plrtrans(node *p, double theta, double lower, double upper)
{
  	int num;
  	double nn, pr, ptheta, angle, angle2, subangle, len;
  	node *pp, *qq;

  	nn = p->width;
  	subangle = (upper - lower) / nn;
  	qq = p;
  	pp = p->next;
  	if (p->tip)
    	return;
  	angle = upper;
  	do 
	{
	    angle -= pp->back->width / 2.0 * subangle;
	    pr = p->r;
	    ptheta = p->theta;
	    if (regular) 
		{
      		num = 1;
      		while (num * subangle < 2 * pi)
        	num *= 2;
      		if (angle >= 0.0)
        		angle2 = 2 * pi / num * (long)(num * angle / (2 * pi) + 0.5);
      		else
        		angle2 = 2 * pi / num * (long)(num * angle / (2 * pi) - 0.5);
   		} 
		else
      		angle2 = angle;
    	if (uselengths)
      		len = fabs(pp->back->oldlen);
    	else
      		len = 1.0;
    	pp->back->r = sqrt(len * len + pr * pr + 2 * len * pr * cos(angle2 - ptheta));
    	if (fabs(pr * cos(ptheta) + len * cos(angle2)) > epsilon)
      		pp->back->theta = atan((pr * sin(ptheta) + len * sin(angle2)) /
                 (pr * cos(ptheta) + len * cos(angle2)));
    	else if (pr * sin(ptheta) + len * sin(angle2) >= 0.0)
      		pp->back->theta = pi / 2;
    	else
      		pp->back->theta = 1.5 * pi;
    	if (pr * cos(ptheta) + len * cos(angle2) < -epsilon)
      		pp->back->theta += pi;
    	if (!pp->back->tip)
      		plrtrans(pp->back, pp->back->theta,
            	angle - pp->back->width * subangle / 2.0,
                angle + pp->back->width * subangle / 2.0);
    	else
      		pp->back->oldtheta = angle2;
    	angle -= pp->back->width / 2.0 * subangle;
    	pp = pp->next;
  	} while (((p != root) && (pp != qq)) || ((p == root) && (pp != p->next)));
}
/** added - danieyek 990208 
 * Find the median of the distance; used to compute the angle to rotate
 * in proportion to the size of the graph and forces.
 * It is assumed that pRootSubNode is also the pivot (during the first 
 * call to this function) - the center, with
 * respect to which node the distances are calculated.
 * If there are only 3 element, element #2 is returned, ie. (2+1)/2.
 * This function now finds the median of distances of all nodes, not only
 * the leafs! 
 */
double medianOfDistance(node *pRootSubNode, boolean firstRecursiveCallP)
{
  	node *pSubNode;
  	double xDelta, yDelta, distance;
  	long i, j;
  	struct dblLinkNode
  	{
    	double value;
    	/* Implement reverse Linked List */
    	struct dblLinkNode *pBack;
  	} *pLink, *pBackElem, *pMidElem, *pFrontElem, junkLink;
  	/* must use static to retain values over calls */
  	static node *pReferenceNode;
  	static long count;
  	static struct dblLinkNode *pFrontOfLinkedList;

  	/* Remember the reference node so that it doesn't have to be passed 
     around in the function parameter. */
  	if (firstRecursiveCallP == true)
  	{
    	pReferenceNode = pRootSubNode;
    	pFrontOfLinkedList = NULL;
   	 	count = 0;
  	}
  	pSubNode = pRootSubNode;
  	/* visit the rest of the branches of current node; the branch attaches to 
     the current subNode may be visited in the code further down below. */
  	while (pSubNode->next != NULL && pSubNode->next != pRootSubNode)
  	{
    	pSubNode = pSubNode->next;
    	if ( pSubNode->back != NULL) 
      		medianOfDistance(pSubNode->back, false);
  	}
  	/* visit this branch; You need to visit it for the first time - at root 
     only! use pRootSubNode instead of pSubNode here because pSubNode stop 
     short just before pRootSubNode (the entry node) */
  	if ( firstRecursiveCallP == true && pRootSubNode->back != NULL)
    	medianOfDistance(pRootSubNode->back, false);
  	/* Why only leafs count?  Modifying it!  */
  	xDelta = nodep[pSubNode->index-1]->xcoord - 
    	nodep[pReferenceNode->index-1]->xcoord;
  	yDelta = nodep[pSubNode->index-1]->ycoord - 
    	nodep[pReferenceNode->index-1]->ycoord;
  	distance = sqrt( xDelta*xDelta + yDelta*yDelta );
  	/* Similar to pushing onto the stack */
  	pLink = (struct dblLinkNode*) Malloc( sizeof(struct dblLinkNode) );
  	if (pLink == NULL)
  	{
    	syslog(LOG_ERR,"Fatal ERROR: drawtree - Insufficient Memory in" 
        " medianOfDistance()!\n");
    	exit(-1);
  	}
  	pLink->value = distance;
  	pLink->pBack = pFrontOfLinkedList;
  	pFrontOfLinkedList = pLink;
  	count = count + 1;
  	if (firstRecursiveCallP == true) 
  	{
    	if (count == 0)
    	{
      		return (double)0;
    	}
    	else if (count == 1)
    	{
      		distance = pFrontOfLinkedList->value;
      		free(pFrontOfLinkedList);
      		return distance;
    	}
    	else if (count == 2)
    	{
      		distance = (pFrontOfLinkedList->value + 
          	pFrontOfLinkedList->pBack->value)/(double)2;
      		free(pFrontOfLinkedList->pBack);
      		free(pFrontOfLinkedList);
      		return distance;
    	}
    	else
    	{
      		junkLink.pBack = pFrontOfLinkedList;
	
      		/* SORT first - use bubble sort; we start with at least 3 elements here. */
      		/* We are matching backward when sorting the list and comparing MidElem and 
         	BackElem along the path;  junkLink is there just to make a symmetric
         	operation at the front end. */
      		for (j = 0; j < count - 1; j++)
      		{
            	pFrontElem = &junkLink;
            	pMidElem = junkLink.pBack;
            	pBackElem = junkLink.pBack->pBack;
        		for (i = j; i < count - 1; i++)
        		{
          			if(pMidElem->value < pBackElem->value)
          			{	
          				/* Swap - carry the smaller value to the root of the linked list. */
            			pMidElem->pBack = pBackElem->pBack;
            			pBackElem->pBack = pMidElem;
            			pFrontElem->pBack = pBackElem;
            			/* Correct the order of pFrontElem, pMidElem, pBackElem and match 
               			one step */
            			pFrontElem = pBackElem;
            			pBackElem = pMidElem->pBack;
          			}
          			else
          			{
            			pFrontElem = pMidElem;
            			pMidElem = pBackElem;
            			pBackElem = pBackElem->pBack;
          			}
        		}
        		pFrontOfLinkedList = junkLink.pBack;
      		}
      		/* Sorted; now get the middle element. */
      		for (i = 1; i < (count + 1)/(long) 2; i++)
      		{
       		 	/* Similar to Poping the stack */
              	pLink = pFrontOfLinkedList;
        		pFrontOfLinkedList = pLink->pBack;
        		free(pLink);
      		}
    		/* Get the return value!! - only the last return value is the valid one. */
      		distance = pFrontOfLinkedList->value;

      		/* Continue from the same i value left off by the previous for loop!  */
      		for (; i <= count; i++)
      		{
        		/* Similar to Poping the stack */
              	pLink = pFrontOfLinkedList;
        		pFrontOfLinkedList = pLink->pBack;
        		free(pLink);
      		}
    	}
  	}
  	return distance;
}
double dotProduct(double Xu, double Yu, double Xv, double Yv)
{
  	return Xu * Xv + Yu * Yv;
} 
/* Calculate angle between 2 vectors; value returned is always between
   0 and pi.  
   Use vCounterClkwiseU() to get the relative position of the vectors.  */ 
double angleBetVectors(double Xu, double Yu, double Xv, double Yv)
{
  	double dotProd, cosTheta, theta, lengthsProd;
  	dotProd = dotProduct(Xu, Yu, Xv, Yv);
  	lengthsProd = sqrt(Xu*Xu+Yu*Yu) * sqrt(Xv*Xv+Yv*Yv);

  	if (lengthsProd < epsilon)
  	{
    	syslog(LOG_ERR,"drawtree - division by zero in angleBetVectors()!\n"
    	"Xu %f Yu %f Xv %f Yv %f\n", Xu, Yu, Xv, Yv);
    	exit(-1);
  	}
  	cosTheta = dotProd / lengthsProd;

  	if (cosTheta > 1) /* cosTheta will only be > 1 or < -1 due to rounding errors */
    	theta = 0; /* cosTheta = 1 */
  	else if (cosTheta < -1)
    	theta = pi; /* cosTheta = -1 */
  	else 
    	theta = acos(cosTheta);
  	return theta;
}
/* As usual, pToSubNode->back is the angle 
   leftLimit is the max angle you can rotate on the left and 
   rightLimit vice versa.
   *pLeftLimit and *pRightLimit must be initialized to 0; without
   initialization, it would introduce bitter bugs into the program;
   they are initialized in this routine.
*/
void leftRightLimits(node *pToSubNode, double *pLeftLimit, 
	double *pRightLimit)
{
  	/* pPivotNode is nodep[pToSubNode->back->index-1], not pPivotSubNode
     which is just pToSubNode->back! */
  	node *pLeftSubNode, *pRightSubNode, *pPivotNode, *pSubNode;
  	double leftLimit, rightLimit, xToNodeVector, yToNodeVector, xLeftVector, 
    	yLeftVector, xRightVector, yRightVector, lengthsProd;

  	*pLeftLimit = 0;
  	*pRightLimit = 0;
  	/* Make an assumption first - guess "pToSubNode->back->next" is the right 
     and the opposite direction is the left! */
  	/* It shouldn't be pivoted at a left, but just checking. */
  	if (pToSubNode->back->tip == true)
  	{
   		/* Logically this should not happen.  But we actually can return pi 
      	as the limit. */
    	syslog(LOG_ERR,"In leftRightLimits() - Pivoted at a leaf! Unable to " 
        	"calculate left and right limit.\n");
    	exit(-1);
  	}
  	else if (pToSubNode->back->next->next == pToSubNode->back)
  	{
    	*pLeftLimit = 0; /* don't pivot where there is no branching */
    	*pRightLimit = 0;
    	return;
  	}
  	/* Else, do this */
  	pPivotNode = nodep[pToSubNode->back->index-1];

  	/* 3 or more branches - the regular case. */
  	/* First, initialize the pRightSubNode - non-repeative portion of the code */
  	pRightSubNode = pToSubNode->back;
  	pLeftSubNode = pToSubNode->back;
  	xToNodeVector = nodep[pToSubNode->index-1]->xcoord - pPivotNode->xcoord;
  	yToNodeVector = nodep[pToSubNode->index-1]->ycoord - pPivotNode->ycoord;

  	/* If both x and y are 0, then the length must be 0; but this check is not
     enough yet, we need to check the product of length also. */
  	if ( fabs(xToNodeVector) < epsilon && fabs(yToNodeVector) < epsilon )
  	{
    	/* If the branch to rotate is too short, don't rotate it. */
    	*pLeftLimit = 0;
    	*pRightLimit = 0;
    	return;
  	}
  	while( nodep[pRightSubNode->index-1]->tip != true )
  	{
    	/* Repeative code */
    	pRightSubNode = pRightSubNode->next->back;

    	xRightVector = nodep[pRightSubNode->index-1]->xcoord - pPivotNode->xcoord;
    	yRightVector = nodep[pRightSubNode->index-1]->ycoord - pPivotNode->ycoord;

    	lengthsProd = sqrt(xToNodeVector*xToNodeVector+yToNodeVector*yToNodeVector)
    	  * sqrt(xRightVector*xRightVector+yRightVector*yRightVector);
   	 	if ( lengthsProd < epsilon )
    	{
      		continue;
    	}
    	rightLimit = angleBetVectors(xToNodeVector, yToNodeVector, xRightVector, 
        	yRightVector);
    	if ( (*pRightLimit) < rightLimit) *pRightLimit = rightLimit;
  	}

  	while( nodep[pLeftSubNode->index-1]->tip != true )
  	{
    	/* First, let pSubNode be 1 subnode after rightSubNode. */
    	pSubNode = pLeftSubNode->next->next;
    	/* Then, loop until the last subNode before getting back to the pivot */
    	while (pSubNode->next != pLeftSubNode)
    	{
      		pSubNode = pSubNode->next;
    	}
    	pLeftSubNode = pSubNode->back; 
    	xLeftVector = nodep[pLeftSubNode->index-1]->xcoord - pPivotNode->xcoord;
    	yLeftVector = nodep[pLeftSubNode->index-1]->ycoord - pPivotNode->ycoord;

    	lengthsProd = sqrt(xToNodeVector*xToNodeVector+yToNodeVector*yToNodeVector)
      	* sqrt(xLeftVector*xLeftVector+yLeftVector*yLeftVector);
    	if ( lengthsProd < epsilon )
    	{
      		continue;
    	}
    	leftLimit = angleBetVectors(xToNodeVector, yToNodeVector, xLeftVector, 
        	yLeftVector);

    	if ( (*pLeftLimit) < leftLimit) *pLeftLimit = leftLimit;
  	}
  	return;
}
/* compute the angle of a vector */
double angleof(double x, double y)
{
  	double theta;

  	if (fabs(x) > epsilon)
    	theta = atan(y / x);
  	else if (y >= 0.0)
      	theta = pi / 2;
    else
      	theta = 1.5 * pi;
  	if (x < -epsilon)
    	theta = pi + theta;
  	while (theta > 2 * pi)
    	theta -= 2 * pi;
  	while (theta < 0.0)
    	theta += 2 * pi;
  	return theta;
}
/* get leftmost and rightmost angle of subtree, put them in node p */
void leftrightangle(node *p, double xx, double yy)
{
  	double firstx, firsty, leftx, lefty, rightx, righty;
  	double langle, rangle;

  	firstx = nodep[p->back->index-1]->xcoord - xx;
  	firsty = nodep[p->back->index-1]->ycoord - yy;
  	leftx = firstx;
  	lefty = firsty;
  	rightx = firstx;
  	righty = firsty;
  	if (p->back != NULL)
    	polartrav(p->back,xx,yy,firstx,firsty,&leftx,&lefty,&rightx,&righty);
  	if ((fabs(leftx) < epsilon) && (fabs(lefty) < epsilon))
    	langle = p->back->oldtheta;
  	else
    	langle = angleof(leftx, lefty);
  	if ((fabs(rightx) < epsilon) && (fabs(righty) < epsilon))
    	rangle = p->back->oldtheta;
  	else
    	rangle = angleof(rightx, righty);
  	while (langle - rangle > 2*pi)
    	langle -= 2 * pi;
  	while (rangle > langle) 
	{
    	if (rangle > 2*pi)
      		rangle -= 2 * pi;
    	else
      		langle += 2 * pi;
  	}
  	while (langle > 2*pi) 
	{
    	rangle -= 2 * pi;
    	langle -= 2 * pi;
  	}
  	p->lefttheta = langle;
  	p->righttheta = rangle;
}
/** 
 * Read comment for totalForceOnNode
 * It supposed to return a positive value to indicate that it has a positive 
 * moment; and negative return value to indicate negative moment. 
 * force perpendicular at norminal distance 1 is taken to be 1.  
 * medianDistance is the median of Distances in this graph. 
 *
                        / Force
                       /  
          |  ToNode   o  > alpha
          |              \  
   yDelta |               \     theta = pi/2 + alpha
          |                  beta = vector (or angle) from Pivot to ToNode
  Pivot   o-----------
             xDelta

  alpha = theta + beta
 */
double forcePerpendicularOnNode(node *pPivotSubNode, node *pToSubNode,
	double medianDistance)
{
  	double totalForce, forceAngle, xDelta, yDelta;
  	double alpha, theta, forcePerpendicular, sinForceAngle, cosForceAngle;

  	totalForce = (double)0;
  	forceAngle = (double)0;
  	totalForceOnNode(pPivotSubNode, pToSubNode, &totalForce, &forceAngle, 
      	medianDistance);
  	xDelta = nodep[pToSubNode->index-1]->xcoord - 
    	nodep[pPivotSubNode->index-1]->xcoord;
  	yDelta = nodep[pToSubNode->index-1]->ycoord - 
    nodep[pPivotSubNode->index-1]->ycoord;
  	/* Try to avoid the case where 2 nodes are on top of each other. */
  	/*
  	if (xDelta < 0) tempx = -xDelta;
  	else tempx = xDelta;
  	if (yDelta < 0) tempy = -yDelta;
  	else tempy = yDelta;
  	if (tempx < epsilon && tempy < epsilon)
 	{
    	return;
  	}
  	*/
  	sinForceAngle = sin(forceAngle);
  	cosForceAngle = cos(forceAngle);
  	theta = angleBetVectors(xDelta, yDelta, cosForceAngle, sinForceAngle);
  	if (theta > pi/2)
  	{
    	alpha = theta - pi/2;
  	}
  	else
  	{
    	alpha = pi/2 - theta;
  	}
  	forcePerpendicular = totalForce * cos(alpha);

  	if (forcePerpendicular < -epsilon)
  	{
    	syslog(LOG_ERR,"drawtree - forcePerpendicular applied at an angle should" 
        " not be less than zero (in forcePerpendicularOnNode()). alpha = %f\n", alpha);
    	exit(-1);
  	}
  	/* correct the sign of the moment */
  	forcePerpendicular = signOfMoment(xDelta, yDelta, cosForceAngle, 
      	sinForceAngle) 
    	* forcePerpendicular;
  	return forcePerpendicular;
}
void improveNodeAngle(node *pToNode, double medianDistance)
{
  	double forcePerpendicular, distance, 
    xDistance, yDistance, angleRotate, sinAngleRotate, cosAngleRotate, 
    norminalDistance, leftLimit, rightLimit, limitFactor;
  	node *pPivot;

  	/* Limit factor determinte how close the rotation can approach the 
     absolute limit before colliding with other branches */
  	limitFactor = (double)4 / (double)5;
  	pPivot = pToNode->back;

  	xDistance = nodep[pPivot->index-1]->xcoord - nodep[pToNode->index-1]->xcoord;
  	yDistance = nodep[pPivot->index-1]->ycoord - nodep[pToNode->index-1]->ycoord;
  	distance = sqrt( xDistance*xDistance + yDistance*yDistance  );

  	/* convert distance to absolute value and test if it is zero */
  	if ( fabs(distance) < epsilon) 
  	{
    	angleRotate = (double)0;
  	}
  	else 
  	{
    	leftRightLimits(pToNode, &leftLimit, &rightLimit);
    	norminalDistance = distance / medianDistance;
    	forcePerpendicular = forcePerpendicularOnNode(pPivot, pToNode,
        	medianDistance);
    	angleRotate = forcePerpendicular / norminalDistance;
    	/* Limiting the angle of rotation */
    	if ( angleRotate > 0 && angleRotate > limitFactor * leftLimit)
    	{
      		/* Left */
      		angleRotate = limitFactor * leftLimit;
    	}
    	else if ( -angleRotate > limitFactor * rightLimit )
    	/* angleRotate < 0 && */ 
    	{
      		/* Right */
      		angleRotate = - limitFactor * rightLimit;
    	}
  	}
  	angleRotate = (double).1 * angleRotate;

  	sinAngleRotate = sin(angleRotate);
  	cosAngleRotate = cos(angleRotate);

  	tilttrav(pToNode,  
    	&(nodep[pPivot->index - 1]->xcoord), 
        &(nodep[pPivot->index - 1]->ycoord), 
        &sinAngleRotate, &cosAngleRotate);

  	polarizeABranch(pToNode,
        &(nodep[pPivot->index - 1]->xcoord),
        &(nodep[pPivot->index - 1]->ycoord));
}
/** 
 * improvtrav for n-body. 
 * POPStack is the stack that is currently being used (popped); 
 * PUSHStack is the stack that is for the use of the next round 
 * (is pushed now) 
 */
void improvtravn(node *pStartingSubNode)
{
	stackElemType *pPUSHStackTop, *pPOPStackTop, *pTempStack;
  	node *pSubNode, *pBackStartNode, *pBackSubNode;
  	double medianDistance;
  	int noOfIteration;

  	/* Stack starts with no element on it */
  	pPUSHStackTop = NULL;
  	pPOPStackTop = NULL;

  	/* Get the median to relate force to angle proportionally. */
  	medianDistance = medianOfDistance(root, true);

  	/* Set max. number of iteration */
  	for ( noOfIteration = (long)0; noOfIteration < maxNumOfIter; noOfIteration++) 
	{
    	/* First, push all subNodes in the root node onto the stack-to-be-used
       	to kick up the process */
    	pSubNode = pStartingSubNode;
    	pushNodeToStack(&pPUSHStackTop, pSubNode);
    	while(pSubNode->next != pStartingSubNode)
    	{
      		pSubNode = pSubNode->next;
      		pushNodeToStack(&pPUSHStackTop, pSubNode);
    	} 
    	while (true)
    	{
      		/* Finishes with the current POPStack; swap the function of the stacks if 
               PUSHStack is not empty */
      		if (pPUSHStackTop == NULL)
      		{
              	/* Exit infinity loop here if empty. */
              	break;
      		}
      		else
      		{
              	/* swap */
              	pTempStack = pPUSHStackTop;
              	pPUSHStackTop = pPOPStackTop;
              	pPOPStackTop = pTempStack;
      		}
      
      		while (pPOPStackTop != NULL)
      		{
           		/* We always push the pivot subNode onto the stack!  That's when we
            	pop that pivot subNode, subNode.back is the node we apply the force
            	to (ToNode). Also, when we pop a pivot subNode, always push all 
            	pivot subNodes in the same ToNode onto the stack. */
              	popNodeFromStack(&pPOPStackTop, &pSubNode);
      
              	pBackStartNode = pSubNode->back;
              	if (pBackStartNode->tip == true)
              	{
                	/* tip indicates if a node is a leaf */
                	improveNodeAngle(pSubNode->back, medianDistance);
      
              	}
              	else
              	{
                	/* Push all subNodes in this pSubNode->back onto the
                	 * stack-to-be-used, after poping a pivot subNode. If
                	 * pSubNode->back is a leaf, no push on stack. */ 
                	pBackSubNode = pBackStartNode;
                	/* Do not push this pBackStartNode onto the stack!  Or the
                	 * process will never stop. */
                	while(pBackSubNode->next != pBackStartNode)
                	{
                  		pBackSubNode = pBackSubNode->next;
                  		pushNodeToStack(&pPOPStackTop, pBackSubNode);
                	} 
                	/* improve the node even if it is not a leaf */
                	improveNodeAngle(pSubNode->back, medianDistance);
              	}
      		}
    	}
  	}
}
/* traverse tree trying different tiltings at each node */
void improvtrav(node *p)
{
  	double xx, yy, cosphi, sinphi;
  	double langle, rangle, sumrot, olddiff;
  	node *pp, *qq, *ppp;;

  	if (p->tip)
    	return;
  	xx = p->xcoord;
  	yy = p->ycoord;
  	pp = p->next;
  	do 
	{
    	leftrightangle(pp, xx, yy);
    	pp = pp->next;
  	} while ((pp != p->next));
  	if (p == root) 
	{
    	pp = p->next;
    	do 
		{
      		qq = pp;
      		pp = pp->next;
    	} while (pp != root);
    	p->righttheta = qq->righttheta;
    	p->lefttheta = p->next->lefttheta;
  	}
  	qq = p;
  	pp = p->next;
  	ppp = p->next->next;
  	do 
	{
    	langle = qq->righttheta - pp->lefttheta;
    	rangle = pp->righttheta - ppp->lefttheta; 
    	while (langle > pi)
      		langle -= 2*pi;
    	while (langle < -pi)
      		langle += 2*pi;
    	while (rangle > pi)
      		rangle -= 2*pi;
    	while (rangle < -pi)
      		rangle += 2*pi;
    	olddiff = fabs(langle-rangle);
    	sumrot = (langle - rangle) /2.0;
    	if (sumrot > langle)
      		sumrot = langle;
    	if (sumrot < -rangle)
      		sumrot = -rangle;
    	cosphi = cos(sumrot);
    	sinphi = sin(sumrot);
    	if (p != root) 
		{
      		if (fabs(sumrot) > maxchange)
        		maxchange = fabs(sumrot);
      		pp->back->oldtheta += sumrot;
      		tilttrav(pp->back,&xx,&yy,&sinphi,&cosphi);
      		polarize(pp->back,&xx,&yy);
      		leftrightangle(pp, xx, yy);
      		langle = qq->righttheta - pp->lefttheta;
      		rangle = pp->righttheta - ppp->lefttheta; 
      		while (langle > pi)
        		langle -= 2*pi;
      		while (langle < -pi)
        		langle += 2*pi;
      		while (rangle > pi)
        		rangle -= 2*pi;
      		while (rangle < -pi)
        		rangle += 2*pi;
      		while ((fabs(langle-rangle) > olddiff) && (fabs(sumrot) > 0.01)) 
			{
        		sumrot = sumrot /2.0;
        		cosphi = cos(-sumrot);
        		sinphi = sin(-sumrot);
        		pp->back->oldtheta -= sumrot;
        		tilttrav(pp->back,&xx,&yy,&sinphi,&cosphi);
        		polarize(pp->back,&xx,&yy);
        		leftrightangle(pp, xx, yy);
        		langle = qq->righttheta - pp->lefttheta;
        		rangle = pp->righttheta - ppp->lefttheta; 
        		if (langle > pi)
          			langle -= 2*pi;
        		if (langle < -pi)
          			langle += 2*pi;
        		if (rangle > pi)
          			rangle -= 2*pi;
        		if (rangle < -pi)
          			rangle += 2*pi;
      		}
    	}
    	qq = pp;
    	pp = pp->next;
    	ppp = ppp->next;
  	} while (((p == root) && (pp != p->next)) || ((p != root) && (pp != p)));
  	pp = p->next;
  	do 
	{
    	improvtrav(pp->back);
    	pp = pp->next;
  	}
  	while (((p == root) && (pp != p->next)) || ((p != root) && (pp != p)));
}
/* get width and depth beyond each node */
void getwidth(node *p)
{
  	double nw, nd;
  	node *pp, *qq;

  	nd = 0.0;
  	if (p->tip)
    	nw = 1.0;
  	else 
	{
	    nw = 0.0;
    	qq = p;
    	pp = p->next;
    	do 
		{
      		getwidth(pp->back);
      		nw += pp->back->width;
      		if (pp->back->depth > nd)
        		nd = pp->back->depth;
      		pp = pp->next;
    	} while (((p != root) && (pp != qq)) || ((p == root) && (pp != p->next)));
  	}
  	p->depth = nd + p->length;
  	p->width = nw;
}
/* compute x and y coordinates */
void coordtrav(node *p, double *xx, double *yy)
{
  	node *pp;

  	if (!p->tip) 
	{
    	pp = p->next;
    	while (pp != p) 
		{
      		coordtrav(pp->back, xx, yy);
      		pp = pp->next;
      		if (p == root)
        		coordtrav(p->back, xx, yy);
    	}
  	}
  	(*xx) = p->r * cos(p->theta);
  	(*yy) = p->r * sin(p->theta);
  	if ((*xx) > maxx)
    	maxx = (*xx);
  	if ((*xx) < minx)
    	minx = (*xx);
  	if ((*yy) > maxy)
    	maxy = (*yy);
  	if ((*yy) < miny)
    	miny = (*yy);
  	p->xcoord = (*xx);
  	p->ycoord = (*yy);
}
/* use angles calculation to improve node coordinate placement */
void coordimprov(double *xx, double *yy)
{
  	int i;

  	if (nbody)       /* n-body algorithm */
    	improvtravn(root);
  	else 
	{          /* equal-daylight algorithm */
	    i = 0;
	    do 
		{
      		maxchange = 0.0;
      		improvtrav(root);
      		i++;
    	} while ((i < MAXITERATIONS) && (maxchange > MINIMUMCHANGE));
  	}
}
/* compute coordinates for tree */
void calculate()
{
  	double xx, yy;
  	int i;
  	double nttot, fontheight, labangle=0, top, bot, rig, lef;

  	for (i = 0; i < nextnode; i++)
    	nodep[i]->width = 1.0;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->xcoord = 0.0;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->ycoord = 0.0;
  	if (!uselengths) 
	{
    	for (i = 0; i < nextnode; i++)
      		nodep[i]->length = 1.0;
  	} 
	else 
	{
    	for (i = 0; i < nextnode; i++)
      	nodep[i]->length = fabs(nodep[i]->oldlen);
 	}
  	getwidth(root);
  	nttot = root->width;
  	for (i = 0; i < nextnode; i++)
    	nodep[i]->width = nodep[i]->width * spp / nttot;
  	if (!improve)
    	plrtrans(root, treeangle, treeangle - ark / 2.0, treeangle + ark / 2.0);
  	else plrtrans(root, treeangle, treeangle - pi, treeangle + pi);
  		maxx = 0.0;
  	minx = 0.0;
  	maxy = 0.0;
  	miny = 0.0;
  	coordtrav(root, &xx,&yy);
  	fontheight = heighttext(fontname);
  	if (labeldirec == fixed)
    	labangle = pi * labelrotation / 180.0;
  	textlength = (double*) Malloc(nextnode*sizeof(double));
  	firstlet = (double*) Malloc(nextnode*sizeof(double));
  	for (i = 0; i < nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		textlength[i] = lengthtext(nodep[i]->nayme, nodep[i]->naymlength,
				fontname);
      		textlength[i] /= fontheight;
      		firstlet[i] = lengthtext(nodep[i]->nayme,1L,fontname)
                      / fontheight;
    	}
  	}
  	if (spp > 1)
    	labelheight = charht * (maxx - minx) / (spp - 1);
  	else
    	labelheight = charht * (maxx - minx);
  	if (improve) 
	{
		coordimprov(&xx,&yy);
		maxx = 0.0;
		minx = 0.0;
		maxy = 0.0;
		miny = 0.0;
		coordtrav(root, &xx,&yy);
  	}
  	topoflabels = 0.0;
  	bottomoflabels = 0.0;
  	rightoflabels = 0.0;
  	leftoflabels = 0.0;
  	for (i = 0; i < nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		if (labeldirec == radial)
        		labangle = nodep[i]->theta;
      		else if (labeldirec == along)
        		labangle = nodep[i]->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		if (cos(labangle) < 0.0 && labeldirec != fixed)
        		labangle -= pi;
      		firstlet[i] = lengthtext(nodep[i]->nayme,1L,fontname)
				/ fontheight;
      		top = (nodep[i]->ycoord - maxy) / labelheight + sin(nodep[i]->oldtheta);
      		rig = (nodep[i]->xcoord - maxx) / labelheight + cos(nodep[i]->oldtheta);
      		bot = (miny - nodep[i]->ycoord) / labelheight - sin(nodep[i]->oldtheta);
      		lef = (minx - nodep[i]->xcoord) / labelheight - cos(nodep[i]->oldtheta);
      		if (cos(labangle) * cos(nodep[i]->oldtheta) +
          		sin(labangle) * sin(nodep[i]->oldtheta) > 0.0) 
			{
        		if (sin(labangle) > 0.0)
          			top += sin(labangle) * textlength[i];
        		top += sin(labangle - 1.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) < 0.0)
          			bot -= sin(labangle) * textlength[i];
        		bot -= sin(labangle - 0.75 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			rig += cos(labangle - 0.75 * pi) * GAP * firstlet[i];
        		else
          			rig += cos(labangle - 1.25 * pi) * GAP * firstlet[i];
        		rig += cos(labangle) * textlength[i];
        		if (sin(labangle) > 0.0)
          			lef -= cos(labangle - 1.25 * pi) * GAP * firstlet[i];
        		else
          			lef -= cos(labangle - 0.75 * pi) * GAP * firstlet[i];
      		} 
			else 
			{
        		if (sin(labangle) < 0.0)
          			top -= sin(labangle) * textlength[i];
        		top += sin(labangle + 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			bot += sin(labangle) * textlength[i];
        		bot -= sin(labangle - 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) > 0.0)
          			rig += cos(labangle - 0.25 * pi) * GAP * firstlet[i];
        		else
          			rig += cos(labangle + 0.25 * pi) * GAP * firstlet[i];
        		if (sin(labangle) < 0.0)
          			rig += cos(labangle) * textlength[i];
        		if (sin(labangle) > 0.0)
          			lef -= cos(labangle + 0.25 * pi) * GAP * firstlet[i];
        		else
          			lef -= cos(labangle - 0.25 * pi) * GAP * firstlet[i];
        		lef += cos(labangle) * textlength[i];
      		}
      		if (top > topoflabels)
        		topoflabels = top;
      		if (bot > bottomoflabels)
        		bottomoflabels = bot;
      		if (rig > rightoflabels)
        		rightoflabels = rig;
      		if (lef > leftoflabels)
        		leftoflabels = lef;
    	}
  	}
  	topoflabels *= labelheight;
  	bottomoflabels *= labelheight;
  	leftoflabels *= labelheight;
  	rightoflabels *= labelheight;
}
/* compute coordinates of tree for plot or preview device */
void rescale()
{
  long i;
  double treeheight, treewidth, extrax, extray, temp;

  treeheight = maxy - miny + topoflabels + bottomoflabels;
  treewidth = maxx - minx + rightoflabels + leftoflabels;
  if (grows == vertical) {
    if (!rescaled)
      expand = bscale;
    else {
      expand = (xsize - 2 * xmargin) / treewidth;
      if ((ysize - 2 * ymargin) / treeheight < expand)
        expand = (ysize - 2 * ymargin) / treeheight;
    }
    extrax = (xsize - 2 * xmargin - treewidth * expand) / 2.0;
    extray = (ysize - 2 * ymargin - treeheight * expand) / 2.0;
  } else {
    if (!rescaled)
      expand = bscale;
    else {
      expand = (ysize - 2 * ymargin) / treewidth;
      if ((xsize - 2 * xmargin) / treeheight < expand)
        expand = (xsize - 2 * xmargin) / treeheight;
    }
    extrax = (xsize - 2 * xmargin - treeheight * expand) / 2.0;
    extray = (ysize - 2 * ymargin - treewidth * expand) / 2.0;
  }
  for (i = 0; i < (nextnode); i++) {
    nodep[i]->xcoord = expand * (nodep[i]->xcoord - minx + leftoflabels);
    nodep[i]->ycoord = expand * (nodep[i]->ycoord - miny + bottomoflabels);
    if (grows == horizontal) {
      temp = nodep[i]->ycoord;
      nodep[i]->ycoord = expand * treewidth - nodep[i]->xcoord;
      nodep[i]->xcoord = temp;
    }
    nodep[i]->xcoord += xmargin + extrax;
    nodep[i]->ycoord += ymargin + extray;
  }
}
void plotlabels( char *fontname )
{
	int i;
	double compr, dx = 0, dy = 0, labangle, sino, coso, cosl, sinl,
		cosv, sinv, vec;
  	boolean right;
  	node *lp;

  	compr = xunitspercm / yunitspercm;
  	if ( penchange == yes )
    	changepen(labelpen);
  	for (i = 0; i<nextnode; i++) 
	{
    	if (nodep[i]->tip) 
		{
      		lp = nodep[i];
      		labangle = labelrotation * pi / 180.0;
      		if (labeldirec == radial)
        		labangle = nodep[i]->theta;
      		else if (labeldirec == along)
        		labangle = nodep[i]->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		if (cos(labangle) < 0.0)
        		labangle -= pi;
			sino = sin(nodep[i]->oldtheta);
			coso = cos(nodep[i]->oldtheta);
			cosl = cos(labangle);
			sinl = sin(labangle);
      		right = ((coso*cosl+sino*sinl) > 0.0) || (labeldirec == middle);
      		vec = sqrt(1.0+firstlet[i]*firstlet[i]);
      		cosv = firstlet[i]/vec;
      		sinv = 1.0/vec;
      		if (labeldirec == middle) 
			{
          		if ((textlength[i]+1.0)*fabs(tan(nodep[i]->oldtheta)) > 2.0) 
				{
            		dx = -0.5 * textlength[i] * labelheight * expand;
            		if (sino > 0.0) 
					{
              			dy = 0.5 * labelheight * expand;
              			if (fabs(nodep[i]->oldtheta - pi/2.0) > 1000.0)
                			dx += labelheight * expand / (2.0*tan(nodep[i]->oldtheta));
            		} 
					else 
					{
              			dy = -1.5 * labelheight * expand;
              			if (fabs(nodep[i]->oldtheta - pi/2.0) > 1000.0)
                			dx += labelheight * expand / (2.0*tan(nodep[i]->oldtheta));
            		}
          		}
          		else 
				{
            		if (coso > 0.0) 
					{
              			dx = 0.5 * labelheight * expand;
              			dy = (-0.5 + (0.5*textlength[i]+0.5)*tan(nodep[i]->oldtheta))
                    		* labelheight * expand;
            		}
            		else 
					{
              			dx = -(textlength[i]+0.5) * labelheight * expand;
              			dy = (-0.5 - (0.5*textlength[i]+0.5)*tan(nodep[i]->oldtheta))
                    		* labelheight * expand;
            		}
          		}
      		} 
			else 
			{
        		if (right) 
				{
					dx = labelheight * expand * coso;
					dy = labelheight * expand * sino;
					dx += labelheight * expand * 0.5 * vec * (-cosl*cosv+sinl*sinv);
					dy += labelheight * expand * 0.5 * vec * (-sinl*cosv-cosl*sinv);
        		} 
				else 
				{
					dx = labelheight * expand * coso;
					dy = labelheight * expand * sino;
					dx += labelheight * expand * 0.5 * vec * (cosl*cosv+sinl*sinv);
					dy += labelheight * expand * 0.5 * vec * (sinl*cosv-cosl*sinv);
					dx -= textlength[i] * labelheight * expand * cosl;
					dy -= textlength[i] * labelheight * expand * sinl;
        		}
      		}
      		plottext( lp->nayme, lp->naymlength,
				labelheight * expand * xscale / compr, compr,
				xscale * (lp->xcoord + dx + xoffset),
				yscale * (lp->ycoord + dy + yoffset), -180 * labangle / pi,
				fontname );
		}
	}
	if (penchange == yes)
		changepen(treepen);
} 
/**
 * Plot part or all of tree on the plotting device 
 */
void plottree( node *p, node *q )
{
	double x1, y1, x2, y2;
	node *pp;

	x2 = xscale * (xoffset + p->xcoord);
	y2 = yscale * (yoffset + p->ycoord);
	if (p != root) 
	{
		x1 = xscale * (xoffset + q->xcoord);
		y1 = yscale * (yoffset + q->ycoord);
		plot(penup, x1, y1);
		plot(pendown, x2, y2);
  	}
  	if (p->tip)
    	return;
  	pp = p->next;
  	do 
	{
		plottree(pp->back, p);
		pp = pp->next;
  	} while (((p == root) && (pp != p->next)) 
		|| ((p != root) && (pp != p)));
}
/** Initialize parameters */
void initialparms()
{
	paperx = 20.04;
	pagex  = 20.04;
	papery = 28.74;
	pagey  = 28.74;
	plotrparms(spp);
	grows = vertical;
	treeangle = pi / 2.0;
	ark = 2 * pi;
	regular = false;
	rescaled = true;
	bscale = 1.0;
	labeldirec = middle;
	xmargin = 0.08 * xsize;
	ymargin = 0.08 * ysize;
	labelrotation = 0.0;
	plotter = DEFPLOTTER;
	hpmargin = 0.02*pagex;
	vpmargin = 0.02*pagey;
	labelavoid = false;
}
/**
 * Initializes a node
 */
void initdrawtreenode( node **p, node **grbg, node *q, int len,
 	int nodei, int *ntips, int *parens, initops whichinit,
	pointarray treenode, pointarray nodep, char *str, char *ch,
	char *inTree, int *i )
{
	int j;
  	boolean minusread;
  	double valyew, divisor;

  	switch ( whichinit ) 
	{
	  	case bottom:
			gnu(grbg, p);
			(*p)->index = nodei;
			(*p)->tip = false;
			for (j=0;j<MAXNCH;j++)
			  (*p)->nayme[j] = '\0';
			nodep[(*p)->index - 1] = (*p);
			break;
	  	case nonbottom:
			gnu(grbg, p);
			(*p)->index = nodei;
			break;
	  	case tip:
			(*ntips)++;
			gnu(grbg, p);
			nodep[(*ntips) - 1] = *p;
			setupnode(*p, *ntips);
			(*p)->tip        = true;
			(*p)->naymlength = len ;
			strncpy ((*p)->nayme, str, MAXNCH);
			break;
	  	case length:
			processlength(&valyew, &divisor, ch,
			  &minusread, inTree, i, parens);
			if (!minusread)
			  (*p)->oldlen = valyew / divisor;
			else
			  (*p)->oldlen = fabs(valyew/divisor);
			if ((*p)->oldlen < epsilon)
			  (*p)->oldlen = epsilon;
			if ((*p)->back != NULL)
			  (*p)->back->oldlen = (*p)->oldlen;
			break;
	  	case hsnolength:
			haslengths = false;
			break;
	  	default:        /* cases hslength,iter,treewt,unitrwt        */
			break;        /* should not occur                        */
  	}
}
/**
 * Prepare for drawing the tree
 */
static void setUpEnvironment()
{
	node *q, *r;
  	boolean firsttree;

	grbg = NULL;
  	pointarray treenode = NULL;   
  	if ( strlen(inTreePath) > 0 )
	{
		intree = fopen( inTreePath, "r" );
	}
	else
	{
		intree = stdin;
	}
  	if ( intree != NULL )
	{
		char *treeContents;
		readTreeContents( intree, &treeContents );
		firsttree = true;
		allocate_nodep( &nodep, treeContents, &spp );
  		treeread( treeContents, &root, treenode, &goteof, &firsttree, nodep, &nextnode,
      		&haslengths, &grbg, initdrawtreenode, true, -1 );
		if ( !haslengths )
			uselengths = false;
		if ( intree != stdin )
			fclose( intree );
		q = root;
  		r = root;
		while (!(q->next == root))
			q = q->next;
  		q->next = root->next;
  		root = q;
  		chuck(&grbg, r);
  		nodep[spp] = q;
  		where = root;
  		rotate = true;
  		initialparms();
	}
	else
		syslog(LOG_ERR,"Couldn't open file %s\n",inTreePath );
}

/**
 * Parse the argument list and set globals
 * Currently: -i for intree (path), -p for plotfile (path)
 * and -n for max iterations (integer), -f for a postscript 
 * font-name No arguments are required. If unspecified, the 
 * infile is stdin and the outfile is stdout. Default font
 * is Times-Roman.
 * @return 1 if the args were OK, 0 otherwise
 */
static int readArgs( int argc, char ** argv )
{
	int i;
	int temp;
	for ( i=1;i<argc;i++ )
	{
		if ( strlen(argv[i]) > 1 )
		{
			if ( argv[i][0] == '-' && argc > i+1 )
			{
				switch ( argv[i][1] )
				{
					case 'i':
						strncpy(inTreePath,argv[i+1],FILE_NAME_LEN);
						break;
					case 'p':
						strncpy(plotFilePath,argv[i+1],FILE_NAME_LEN);
						break;
					case 'n':
						sscanf( argv[i+1], "%i", &maxNumOfIter );
						if ( maxNumOfIter > 0 )
							maxNumOfIter = 50;
						break;
					case 'f':
						strncpy( fontname, argv[i+1], FONT_NAME_LEN );
						break;
					case 'B':
						if (strcmp(argv[i+1],"true")==0)
      						uselengths = true;
    					else
							uselengths = false;
						break;
					case 'C':
						sscanf(argv[i+1],"%lf%*[^\n]", &charht);
						break;
					case 'I':
						sscanf(argv[i+1],"%d", &temp);
						// no improvement
						if ( temp == 0 )
						{
							nbody = false;
							improve = false;
						}
						// equal daylight
						else if ( temp == 1 )
						{
							nbody = false;
							improve = true;
						}
						// nbody
						else if ( temp == 2 )
						{
							nbody = true;
							improve = true;
						}
						break;
					default:
						return 0;
						break;
				}
			}
		}
	}
	return 1;
}
/**
 * Set sensible defaults for everything
 */
static void setDefaults()
{
	inTreePath[0] = 0;
	plotFilePath[0] = 0;
	maxNumOfIter = 50;
	strncpy( fontname, "Times-Roman", FONT_NAME_LEN );
	uselengths = haslengths;
	charht = 0.3333;
	improve = true;
	nbody = false;
}
/** 
 * Originally: loop to make preview window and decide what to do with it 
 * Intially cut out but then realised that it made initialisations.
 */
void user_loop()
{
    xscale = xunitspercm;
    yscale = yunitspercm;
    plotrparms(spp);
    calculate();
    rescale();
}
/**
 * Main entry point. Completely rewritten because we don't want 
 * (a) dependencies on non-ANSI-C software such as X11, and (b) 
 * to remove the user interaction. Just a plain old commandline 
 * interface with options works best for integration with an 
 * external GUI. Also we will read from stdin and write to stdout 
 * so we can pipe this command in and out of others such as 
 * ghostscript without writing any files.
 */
int main( int argc, char **argv )
{
	setDefaults();
	if ( readArgs(argc,argv) )
	{
		setUpEnvironment();
		user_loop();
		if ( strlen(plotFilePath)>0 )
		{
			plotfile = fopen( plotFilePath, "w" );
			if ( plotfile == NULL )
			{
				syslog(LOG_ERR,"Couldn't open %s for writing\n", 
					plotFilePath );
				plotfile = stdout;
			}
		}
		else
		{
			plotfile = stdout;
		}
		initplotter();
		numlines = 1;
		drawit( fontname, &xoffset, &yoffset, root );
    	finishplotter();
    	fclose( plotfile );
	}
	else
		syslog(LOG_ERR,"Drawtree: invalid argument list\n");
	return 0;
}
